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Abstract 


Two mitogenomes of Cladiucha were newly reported and showed typical pattern of gene arrangement. 
The phylogenetic position of Cladiucha was obtained from tree reconstruction using various data treatment 
methods and substitution models. Significant heterogeneity in the nucleotide composition and mutational 
biases was found in the mitochondrial protein-coding genes, and the third codon position exhibited high levels 
of saturation. Therefore, 14 datasets were conducted under both site-homogeneous and site-heterogeneous 
models. The following conclusions were drawn from the phylogenetic analyses: (i) the monophyly of 
Tenthredinidae was confirmed, (ii) the monophyly of Allantinae + Tenthredininae + Megabelesinae was 
approved, and (iii) within the family, ((((Tenthredininae + Fenusinae) + Allantinae) + Megabelesinae) + 
Nematinae) is probably the most acceptable cladogram for the phylogeny of Tenthredinidae, which is also 
supported by the morphological analysis and a comparative study on the rRNA secondary structure. 
Divergence time estimation analyses indicated that diversification of the major superfamilies of the suborder 
Symphyta occurred around 232.9 Ma, and the splits of Tenthredinidae were dated to 146 Ma, which 
corresponds to the origin of the earliest lineages of flowering plants and major diversifications of core 
angiosperms, respectively. The Cladiucha arose in the Mid-Miocene; at that time, magnoliids are rapidly 


undergoing genus-species differentiation. 
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1. Introduction 


Tenthredinidae is the largest and a considerably complex family of the paraphyletic suborder 
Symphyta (Hymenoptera) with 5645 extant species in 400 general". This family has been variously divided 
into five, six, seven), eight'®7!, 103 and 11"! subfamilies, or into four families and 12 subfamilies H% 
or six families and 17 subfamilies!!!!, The numbers of the subfamilies within Tenthredinidae and their 
contents are not quite consistent, though some tribes have been recognized in most systematic research 
studies. Megabeleses Takeuchi, 1952! and Cladiucha Konow, 1902121 are the only two genera among 
Tenthredinidae members of which the species are known to feeding on plants of Magnoliaceae, a basal 
lineage of Angiosperms. Benson"! set a monotypic tribe, Cladiuchini, for Cladiucha Konow under 
Emphytinae (=Allantinae). Takeuchi“! and Abe & Smith" placed the two genera into Allantinae of 
Tenthredinidae. WeiU?! pointed out that Megabeleses, Cladiucha, and two additional new genera 
represented a unique lineage among the Tenthredinidae and set a subfamily Megabelesesinae (= 
Megabelesinae) to place them, which might locate between Allantinae and Tenthredininae. This opinion 
was followed in the new system provided by Wei & Nie", However, Taeger et al. F still placed the four 
genera into Allantinae. This systematic inconsistency indicates the requirement of more detailed studies to 
better understand their evolutionary history by using comprehensive molecular data under the current 
approaches. 

The mitogenome data has recently been a widely applied tool in resolving longstanding questions 
about the evolutionary history of organisms due to their rapidly increasing number conjoint with their 
relative compactness in size and structure, inheritance type, absent or very infrequent recombination, and 
high mutation rate“! 14, But then, the presence of high A+T content, lineage-specific compositional 
heterogeneity, within-site rate variation (heterotachy) as well as nucleotide substitution saturation can lead 
to challenges in phylogeny reconstructions!5-3!, These undesirable features can frequently cause the 
forming of analytical artifacts such as long-branch attraction and deceptive phylogenetic relationships. The 
occurrence of analytical objects is a commonly encountered problem in the phylogeny reconstruction of 
intra- and interlineages of insect orders, most notably in Hymenopteral?!. 

The effects of these biases in tree reconstruction are reduced using a variety of approaches, such as a 
series of data coding schemes and more sophisticated evolutionary model settings. Comparing the results of 
different data coding regimes is widely used to improve the signal-to-noise ratio in mitogenome data". In 
general, misleading signals from protein-coding genes (PCGs) are not produced equally in the first, second, 
and third codon positions. In fact, the evolutionary rate of the third codon position is higher than that in the 
first and second codon positions. In some studies, the third codon positions have been removed from 
datasets. In addition, purine-pyrimidine (RY) coding (A&G—R; C&T— Y) in the third codon position is a 
common data coding to alleviate the effects of saturated synonymous nucleotide substitutions and rate 
heterogeneity among codon positions">!*°l, However, this position may provide useful phylogenetic 


signals in some groups. Removing or replacing sites may confound the phylogenetic relationship to some 


extent!!®!70] Regier et al." and Zwick et al. proposed and designed a new data treatment, Degen 
coding, to reduce nucleotide compositional heterogeneity and improve the resolution of deep-level 
arthropod relationships. This data scheme can retain nonsynonymous changes at the third codon position 
while eliminating all of the synonymous changes by extending other coding schemes to fully degenerate all 
of the codons!! "201, 

In contrast to the wide application of PCGs in constructing phylogenetic relationships, rRNAs are used 
with caution. It has become a consensus, that rRNA could not be the protagonist in the inference, due to the 
ineffectiveness of aligning. Even using the alignments guidance by the secondary structure, few works 
could obtain satisfactory results. Due to insufficient data accumulation, homology structural motif have 
been explored in only a few taxa (21. 

Here, we sequenced and reported two mitogenome of Cladiucha. We compiled 14 concatenate 
datasets and applied several commonly used phylogenetic inference methods based on both heterogeneity 
and homogeneity models to overcome the systematic bias arising from nonstationarity. We also estimated a 
time-frame of Symphyta evolution, especially the origin and diversification of Tenthredinidae. Finally, we 
discussed the application of the autapomorphy of rRNAs in the phylogeny frame. Our objective was to 
investigate the effect of compositional heterogeneity on phylogenetic reconstruction and usage of 


autapomorphy of the secondary structure of rRNAs within the rapid radiation of Tenthredinidae. 


2. Methodology 


DNA library construction and sequencing 

Total DNA was extracted from Cladiucha magnoliae, and C. punctata using an E.Z.N.A.® Tissue 
DNA Kit (Omega, Norcross, GA) and was stored at -20?C according to the manufacturer's 
instructions. Sequencing libraries with approximately 250-bp insertions were constructed using a 
NEXT flex™ Rapid DNA-Seq Kit (Illumina, San Diego, CA) in accordance with the manufacturer's 
protocol. Each library was sequenced using an Illumina HiSeq 4000 to generate 150-bp paired-end 
reads at Shanghai Majorbio Bio-pharm Technology Co., Ltd. The genomic DNA yielded a total of 
1.32G raw reads (SRR9998491) and 6.24G raw reads (SRR11177465). 

Mitochondrial genome assembly and annotation 

The C. magnoliae and С. punctata reads were imported into GENEIOUS ΚΙΙ 


(http://www.geneious.com)4! and assembled into contigs. The mitogenomes of Megabeleses 


magnoliae (unpublished), M. liriodendrovorax (unpublished), Tenthredo tienmushana (KR703581) 
and Allantus luctifer (KJ713152) were used as references with ‘medium-low sensitivity’ parameters. 
M. magnoliae was the primary reference among these mitogenomes when concerning rrnS. MITOS 


(http://mitos.bioinf.uni-leipzig.de/index.py) was used to predict the protein-coding, transfer RNA, and 


ribosomal RNA genes with annotation from a reference mitogenome. The total tRNA genes were 
identified by MITOS?5! using the invertebrate mitochondrial genetic code. Geneious v11.0.3 


(http://www.geneious.com) was used to determine the initiation and termination codons of the PCGs 


referring to the sequences of other symphytan species with subsequent manual adjustment. Both the 3’ 
and 5’ ends of rrnS were annotated by secondary structure comparison, rather than by the flanking 
genes. 

Saturation, nucleotide compositions and heterogeneity 

MEGA v7.0 26 was used to calculate the A+T content of the nucleotide sequences and the 
relative synonymous codon usage (RSCU). In the strand encoding the majority of the PCGs, the strand 
asymmetry was calculated using the following formulas: GC-skew = (G-C) / (G+C) and AT-skew = 
(А-Т) / (A+T). Bubble charts (Fig. 1) were generated using JMP. Three-dimensional scatter plots of 
the A+T- and GC-skews and A+T% of four datasets (Fig. 7) were plotted using SigmaPlot 149. 

The disparity index (Ip) detects differences in the evolutionary patterns between a pair of 
sequences, thereby indirectly measuring the level of base compositional heterogeneity. P values [25] 
smaller than 0.05 were considered significant. The /p values were calculated for six nucleotides 
datasets (PCG123, PCGI2, PCGI, PCG2, and PCG3) and an amino acid dataset (AA) using 1000 
Monte Carlo replications, as implemented in MEGA 7.0 61. 

AliGROOVEP?! was used to analyze the heterogeneity of sequence divergence within various 
concatenated datasets for a default sliding window size. Indels in the nucleotide dataset were treated 
as ambiguities, and the BLOSUM62 matrix was used as the default amino acid substitution matrix. 

The saturation level in the nucleotide and amino acid sequences was visualized and assessed using 
DAMBE ?! and ASaturA P! respectively. 

Secondary structure prediction 

Secondary structures of ribosomal RNA (rrnS and rrnL) were predicted and inferred from 
alignments with models Tenthredo tienmushana (GenBank accession No. KR703581) and Allantus 
luctifer (GenBank accession No. KJ713152). First, the primary sequence and the secondary structure 
of the model species were aligned in MARNA'?! to output a consensus sequence and a consensus 
structure. This output file was then imported into SSU-ALIGN"?! to predict the rRNA secondary 
structures of C. magnoliae and C. punctata. Lastly, minor changes were made to transform the 
primary results into relative secondary structures manually. VARNA v3-93?^! and RnaViz 2.0.30?! 
were used to draw the secondary structures of rrnL and rrnS. Helix numbering, with minor 
modifications, was performed, following the Apis mellifera rRNA secondary structures. the RNA 
structural logo shown in Figures 3 and 4 were generated by RNALogo""l. 

Data sets 

The ingroup taxa for the phylogenetic analyses comprised 40 symphytan species, including eight 
taxa representatives of Tenthredinidea (Supplementary Table 1). Four taxa were also included as 


outgroups. 


A total of 13 PCGs were aligned individually, excluding the stop codons by MAFFT v7P5l in 
PhyloSuite v1.1.5. Amino acid alignment of the 13 PCGs was performed using Clustal X. The two 
ribosomal RNAs were aligned by MAFFT v7. 

To eliminate the effect of saturation and compositional heterogeneity on the phylogenetic 
reconstruction, the following 14 concatenate datasets were compiled: 1) PCG: 13 protein-coding 
genes; 2) PCG12: 13 protein-coding genes excluding third codon positions; 3) PCG12RY: 13 
protein-coding genes with RY coding strategy; 4) PCGDegen: 13 protein-coding genes with Degen 
coding strategy; 5) 9PCG: nine protein-coding genes excluding four saturated genes; 6) 
9PCG+4PCG12RY: nine protein-coding genes and four saturated genes with RY coding strategy; 
and all six data sets above combined with two RNAs: 7) PCG+RNA; 8) PCG12+RNA; 9) 
PCG12RY+RNA; 10) PCGDegen+RNA; 11) 9PCG+RNA; 12) 9PCG+4PCG12RY+RNA; and 
two amino acid data sets: 13) AA: 13 protein-coding genes translated into amino acids; 14) 11AA: 
unsaturated 11 amino acids. 

Phylogenetic analysis 

The phylogenies were estimated using the Maximum Likelihood (ML) and Bayesian Inference 
(BI) methods. The ML analyses were conducted using an IQ-TREE web server 
(http://iqtree.cibiv.univie.ac.at/) and RAxML v8.2.10!°!. The Bayesian analyses were performed using 
MrBayes v3.2.6"!!! and PhyloBayes v3.37), 


Before performing the model-based phylogenetic analyses, we determined the best partitioning 
scheme and model using PartitionFinder ν].1.1145:4] The branch lengths and search algorithm settings 
were "unlinked" and "greedy". The Bayesian information criterion (BIC) scheme and the corrected 
Akaike information criterion (AIC) were selected for the BI and ML methods, respectively. 
Accordingly, the datasets were partitioned into up to 16 subsets, and the best models were used for 
subsequent phylogenetic analyses. 

The ML tree was calculated with branch support estimated from 500 bootstrap replicates. The BI 
analyses consisted of four separate runs of four chains each, for 10 million generations, sampling 
every 1000 trees and discarding the first 25% generations as burn-in. An additional Bayesian analysis 
was undertaken in PhyloBayes 3.3?! under the heterogeneous model CAT-GTR for 9PCG and 
CAT-MtArt for 11AA, including four independent chains. These analyses were terminated after the 
likelihood of the sampled trees had stabilized, and the chains had satisfactorily converged (maxdiff « 
0.3). The phylogenetic trees were visualised by FigTree v1.4.2 
(http://tree.bio.ed.ac.uk/software/figtree/). 


Time estimation 
Bayesian estimation of divergence times was performed using the dataset of 11AA in MCMCTree 
implemented in ће РАМІ. package v4.9"! using the approximate likelihood method"6, BASEML 


was used for the estimation of substitution rate per site, and this rate was used as the prior for 


Bayesian analysis. MCMC was run by 50 x 10,000 iterations with the mtArt substitution model. The 
soft bound of the split of the superfamily Tenthredinoidea from others [between 260-290 million 
years ago (Ma)] was used as external secondary calibration based on previous studies^^4?!, The 
inferred nodal age of holometabolous insects (between 300-360 Ma) was used for the calibration of 
the root 2501, 


All of the related files have been uploaded to Figshare 
(https://figshare.com/account/home#/projects/78339) 


3. Results 


Cladiucha mitogenomes 

The nearly complete mitogenomes of C. magnoliae and C. punctate are 15,214 bp and 15,301 bp 
in length, respectively, and are accessed in GenBank (accession number: MT295305-MT295306). The 
organizations of these two mitogenomes are summarized in Table 1. 

Compared with the putative ancestral mitogenome of insects, only tRNA rearrangement events 
were detected. The position between trnM and trnI was swapped, and trnR was translocated from the 
second position in the trnA-trnR-trnN-trnSI-trnE-trnF cluster to the penultimate position. The 
A+T-rich region was partly sequenced (Fig. 5). 

The gene overlaps, and gene gaps between C. magnoliae and C. punctate occurred in the same 
position. We identified six overlapping positions and 18 intergenic spacers. Twelve pairs of genes 
were directly adjacent. The total overlap length was 16 bp in both C. magnoliae and C. punctata, with 
each overlap ranging from 1 to 7 bp. The longest overlapping nucleotide was found between atp6 and 
atp8. There were a total of 370 bp and 442 bp intergenic spacer sequences in C. magnoliae and C. 
punctata, respectively, which ranged from 1 to 151 bp in size. The longest intergenic spacer sequences 
for C. magnoliae and C. punctata of 130 bp and 151 bp, respectively, were located between trnV and 
rrnS. Using homologous searches, no significant similarities were found among the longest noncoding 
region in C. magnoliae, C. punctata, or other identified symphytan species. 

Nine PCGs are located on the majority strand (J-strand), and four PCGs are located on the 
minority strand (N-strand) among the 13 PCGs of Cladiucha mitogenomes. The total lengths of the 
PCGs of these two species were 11,215 bp and 11,168 bp, respectively, accounting for 73.7% and 
73.0%, respectively, of the entire genomes. All of the PCGs started with ATN codons, including four 
genes (nad2, atp6, cox3, and cob) starting with ATG, five genes (сох/, сох2, αἴρδ, nad5, and nad4L) 
starting with ATT and four genes (nad3, nad4, паб, and nad1) starting with ATA. The stop codons 
were mostly TAA, except for nad4, which had T as the stop codon. 


The nearly complete mitogenomes of both C. magnoliae and C. punctata contain 22 tRNA genes. 
There are 14 and eight tRNAs encoded by the J and N strands, respectively. The position and 
orientation of the predicted tRNAs are identical in the two species, ranging in size from 64 bp (trnD) 
to 75 bp (trnC) (Table 1). Only trnS/ (AGN) was missing a dihydrouridine (DHU) arm in both species. 
The predicted anticodons are identical to those of other known symphytan mitochondrial genomes 
(Fig. 6). 

Both of the rRNA genes were encoded on the N-strand. The rRNAs rrnL were flanked by trnL/ 
and trnV. rrnS were downstream of the trnV, with long intergenic space sequences. The length of the 
rrnL gene was 1344 bp in C. magnoliae and 1341 bp in C. punctata, with the same A+T contents of 
84.7%. The lengths of the rrnS genes for C. magnoliae and C. punctata were 813 bp (84.0% A+T 
content) and 835 bp (84.4% A+T content), respectively. There are four domains with 27 helices in 
rrnS and six domains with 46 helices in rrnL in both species (Fig. 3-4). 

Saturation, compositional heterogeneity and codon usage bias within nt data set 

The nucleotide composition in the mitogenomes was described using two statistics: the A+T 
content and the AT- and GC-skewsP". The AT-skew values of C. magnoliae and C. punctata were 
positive (0.0169 and 0.0145, respectively), whereas the corresponding GC-skew values were negative 
(-0.0183 and -0.1860) for the entire mitogenome. The total A+T content of C. magnoliae was 82.8%, 
ranging from 74.5% (cox1) to 90.1% (atp8), whereas that of C. punctata was 82.8%, ranging from 
74.8% (cox) to 91.4% (atp8). The highest A+T contents of 96.4% and 96.0% were found at the third 
codon position for C. magnoliae and C. punctata, respectively. 

The RSCU of C. magnoliae and C. punctata are presented. A significant correlation has been 
found between codon usage and nucleotide composition in C. magnoliae, C. punctata, and also in 
other insect mitogenomes. UUU-Phe, UUA-Leu, AUU-Ile, AUA-Met, and AAU-Asn had the highest 
usage rate among all amino acids. The highest RSCUs were found for UUA-Leu (5.06 and 5.26 for C. 
magnoliae and C. punctata, respectively). CUC-Leu, CUG-Leu, GUC-Val, UCG-Ser, CCG-Thr, 
ACG-Thr, GCG-Ala, CGC-Arg, CGG-Arg, and AGC-Ser all have a high CG content and were least 
frequently used. Crozier and Crozier'??! discussed how the effect of amino acid occurrence depends on 
the nucleotide composition of codon usage, which can be calculated in terms of the ratio of the total 
occurrence of C+G (Pro, Ala, Arg, and Gly) to A+T (Phe, Ile, Met, Tyr, Asn, and Lys) in the 
mitogenomes. This ratio was 0.18 and 0.25 in C. magnoliae and C. punctata, respectively, which is 
close to that of other symphytan species (0.27-0.3 153541, 

The total A+T% of PCGs + rRNA of all of the included hymenopteran species ranged between 
83.4096 (Taeniogonalos taihorina) and 75.1396 (Xyela curva), with a mean of 79.16 (+2.18)%. 
Among Tenthredinidae, Monocellicampa pruni had the lowest A+T% (76.23%), whereas Cladiucha 
had the highest A+T% (82.25% and 82.15%) (Fig. 2A). Irrespective of whether the third codon 


positions (Fig. 2B) or the four oversaturated PCGs (Fig. 2C) were removed, X. curva and T. taihorina 


had the lowest and highest A+T contents, respectively, in Hymenoptera, and M. pruni, and the two 
Cladiucha species had the lowest and highest A+T contents, respectively, in Tenthredinidae. Counting 
only the 13 PCGs, the hymenopteran species with the lowest A+T content was Orussus occidentalis 
(74.12%), followed by X. curva, with 74.41% A+T (Fig. 2D); otherwise, the same data were obtained 
as for the three abovementioned datasets. 

However, it was difficult to characterize each family based on the A+T content alone. Therefore, 
we drew three-dimensional scatter plots of the content, the AT-skew, and the GC-skew for 44 species 
in the PCG+RNA dataset. A radius was assigned to each coordinate point and used to simulate the 
locations of the relatively close species (Fig. 5). The different taxa occupied relatively independent 
positions in space, especially Xyelidae (Fig. 5B) and Pamphiloidae (Fig. 5C). In Tenthredinidae, 
Cladiucha occupied a relatively independent position in contrast to the other Tenthredinidae species, 
which can be confirmed in Fig. 7. As shown in the left of Fig. 7, compared to the rest Tenthredinidae 
species, Cladiucha had the highest A+T content overall and for each codon. When compared within 
Symphyta, Cladiucha had the highest AT-3%. 

The Ip test showed different evolutionary patterns for all six datasets; that is, evolution was 
significantly nonhomogeneous. 

A total of 946 pairwise comparisons were made for each dataset. The null hypothesis was rejected 
in 800 comparisons at the 5% significance level by comparing all 13 PCGs. 

The number of rejected comparisons based on individual codon partitions suggested that the level 
of base compositional heterogeneity was the lowest in PCG2 (485), followed by PCGI2 (677), PCG1 
(681) and PCG3 (703). The maximum Ip was 19.59 in PCGI (between Taeniogonalos taihorina and 
Anopheles gambiae) and 30.88 in PCG3 (between C. magnoliae and Trachelus iudaicus), suggesting 
higher evolutionary variability for the third codon than the first codon. The maximum Ip of 7.53 for 
the AA dataset was found between 7. taihorina and Anopheles gambiae, which was below that for the 
13 PCGs (27.076, between Paroster microsturtensis and T. taihorina) and confirmed that the amino 
acid characters were less compositionally heterogeneous than nucleotide characters. The most 
homogeneous data set corresponded to PCG2, which had the lowest maximum Ip (6.79, between 
Parapolybia crocea and Neopanorpa phlchra) and the fewest rejected comparisons (485). However, 
we were not able to use the PCG2-based phylogeny results to cover the Pamphiloidae. 

The PCG3 mean ID was higher than that of PCGI in each hymenopteran species, except for three 
Xyelidae species and two apocritan species. (Fig. 7). In addition to having the lowest and highest A+T 
contents, Tenthredinidae, Cladiucha and Monocellicampa had much higher nt123Sum and nt3Sum 
than the other species. Two Xyela species and T. taihorina in Apocrita had both nt12Sum and nt1 Sum 


above 200 as well as the lowest and highest A+T contents in Symphyta, respectively. 


The divergence between the nucleotide datasets (PCG123, PCGI23RNA, PCGI2, PCGI2RNA, 
9PCG123 and 9PCG123RNA) and the amino acid datasets (AA and 11AA) for the 44 species were 
analyzed using AliGROOVE software. 

In the PCG dataset (Fig. 8A), the species with high heterogeneity were mainly distributed in 
Orussidae, Siricidae, and Cephidae, with one species each in Pergidae and Apocrita. However, the 
heterogeneity was reduced to varying degrees in the coding dataset (the PCG datasets with data 
recoding; Fig. 8 B-E). Relatively positive scores were obtained for the amino acid datasets for nearly 
all of the family comparisons. The 11AA datasets were the least likely to violate the phylogenetic 
assumptions. However, with the addition of RNA, except for the taxa above, almost all of the 
Cephidae showed heterogeneity (Fig. 8 F-H). 

Comparative analysis of rRNA 

In the rrnL secondary structure, H837 was a long variable stem with a small loop, whereas H991 
and H1196 had variable helical lengths and loop sizes. The predicted structures of H563, H579, H777, 
H822, H1830, H2023, H2043, H2455, and H2547 were well conserved in C. magnoliae, C. punctata, 
and other reported insect species. H976 was redundant in two species and T. tienmushana. H1775 and 
H2347 were conserved with three pairs and a small loop. In the rrnL secondary structure of C. 
magnoliae and C. punctata, positions 293 and 294 (H777 helices) harbored AA, whereas 7. 
tienmushana and Birmella discoidalisa (GenBank accession number KR703581, MF197548) harbored 
GU and A. luctifer (NC024664) harbored GA at the same site. The U at position 411 in H946 was 
replaced by A; the A-U (position 580 in H579) was changed to U-A; the G (position 1208 in H2588) 
was replaced by A; and the UU (position 1270-1271 in H2646) was changed to CC in T. tienmushana, 
A. luctifer and B. discoidalisa. 

In the rrnS secondary structure, domain III was the most conserved domain in Tenthredinidae. 
H47 was a commonly observed variable loop, whereas H500, H769, H944, H1047, H1399 and H1506 
were highly conserved in the two Cladiucha species and other symphytan species. In terms of 
sequence and structure, H9 and H367 were the most conserved helices in the two Cladiucha species 
and in Т. tienmushana, A. luctifer, and B. discoidalisa. The H960, H577, and H1241 helices in the two 
Cladiucha species were more similar to T. tienmushana than A. luctifer and B. discoidalisa. 

Phylogenetic relationships 

The following conclusions are drawn from the results of 32 schemes (Fig. 9): (1) Xyelidae is the 
root of Hymenoptera, with Macroxyelinae and Xyelinae forming the monophylum. (ii) Within 
Tenthredinoidea, Argidae and Pergidae form a monophylum as sisters to the remaining Cimbicidae 
and Tenthredinidae. (11) Within Tenthredinidae, Nematinae (Analcellicampa xanthosoma + 
Monocellicampa pruni) is a sister group of (Allantinae + Tenthredininae + Megabelesinae + 


Fenusinae) and the basal branch of Tenthredinidae. 


However, the results of the 32 phylogenetic trees are inconclusive on the following issues: (i) the 
monophyly of Pamphilioidea; (ii) the relationships within Unicalcarida; and (iii) although most 
branches are highly supported, seven topologies have emerged in Tenthredinidae. 

Eight schemes supported topology No. 3., which was considered to be optimal in this work and 
was supported by phylogenetic tree construction with only eight Tenthredinidae species. Within 
Tenthredinidae, all 13 PCGs were unsaturated. Therefore, eight taxa_PCG and eight taxa_PCG + 
rRNA datasets were used to build the BI tree, both of which supported topology No. 3. Moreover, this 
topology was confirmed by the rRNA secondary structure. Topology No. 1 confirmed to topology No. 
3, except for the relative position of Megabelesinae and Allantinae. 

Three schemes supported topology No. 5. and three other schemes, including the BI and ML trees 
of 11AA, supported topology No. 6. They both cover the monophyly of Allantinae + Tenthredininae + 
Megabelesinae. However, the secondary structure (Fig. 10) shows that Allantinae, Tenthredininae, and 
Fenusinae have multiple common derivatives, excluding topologies No. 5 and No. 6 to some extent. 

Although there are nine schemes supporting topology No. 4 and six schemes supporting topology 
No. 2, respectively, no apomorphic morphological character was found to confirm this result. 

Dated phylogeny 

A chronogram of divergence times for the included hymenopteran species based on the obtained 
tree topology is shown in Fig. 11. According to the divergence time analysis, the crown age of 
Hymenoptera was estimated as 318.1 Ma [95% CI (300.88 — 349.61 Ma)] corresponding to 
Pennsylvanian (Carboniferous). Diversification of the major superfamilies of Symphyta occurred 
between the beginning of the Late Triassic (Carnian, 232.9 Ma) and Jurassic-Cretaceous transition 
(145.0 Ma). The splits of the families of the superfamily Tenthredinoidea were estimated as 
corresponding to Late Jurassic [Pergidae + Argidae, 159.5 Ma, 95% CI (94.73 — 224.25 Ma)] and 
Early Jurassic [Cimbicidae + Tenthredinidae, 187.5 Ma, 95% CI (137.98 — 238.00 Ma)]. The splits of 
the subfamilies Nematinae and Allantinae from other subfamilies of the family Tenthredinidae were 
dated to 146 Ma [95% CI (102.11 — 195.15 Ma)] and 116.3 Ma [95% CI (77.53 — 162,29 Ma)], 
coinciding to Late Jurassic (Tithonian) and Early Cretaceous (Aptian), respectively. The split time of 
the subfamily Tenthredininae from Heterarthrinae + Megabelesesinae was inferred as 103.6 Ma [95% 
CI (64.01 — 148.06 Ma)], corresponding to Early Cretaceous (Aptian). The divergence between the 
subfamilies Heterarthrinae and Megabelesesinae was estimated to take place in Late Cretaceous 
[Coniacian, 86.4 Ma]. The divergence time of C. magnoliae and C. punctata corresponds to the 


Mid-Miocene [15.5 Ma, 95% (5.7 — 34.39 Ma)]. 


4. Discussions 


Mitogenome characteristics of Cladiucha 

The translocation of trnR is a shared derived character of Cladiucha, which is a novel event in 
Symphyta and could be identified as a synapomorphy of Cladiucha. The trnI, trnQ, and trnM gene 
cluster rearrangement events are too random to reconstruct the pattern of genome rearrangements. 
Thus, the evidence from gene rearrangements could not be used to determine the phylogeny of 
Tenthredinidae. 

In C. magnoliae, C. punctata, a conserved overlap and two conserved noncoding sequences were 
found. The conserved overlap (A) located between αἱρὃ and cox3, and the conserved noncoding 
sequences (AA and T) located between nad4L and trnT, and between trnP, and nad6. These 
characteristics were also found in A. luctifer, T. tienmushana and B. discoidalisa. 

We found 21 and 20 mismatched pairs in the tRNAs of C. magnoliae and C. punctata, 
respectively, that were composed of G-U (10 pairs in C. magnoliae and 11 pairs in C. punctata), U-U 
(9 pairs in C. magnoliae and seven pairs in C. punctata) and G-A (two pairs in each species) 
mispairing, which is common in Hymenoptera!^555?l, The mismatched pairs were mainly located on 
the DHU and the anticodon stems. 

The /p corresponds to the disparity between the observed compositional difference for two 
sequences and the expected difference under homogeneity. That is, the /p is zero when the 
homogeneity assumption is satisfied", 

Homogeneity (/p-0) was observed between Cladiucha and Tremex, Trichiosoma, and 
Megalodontes (the first and second datasets). However, these taxa seldom share biological or 
evolutionary characteristics. Cladiucha distributes in southern China and neighboring countries of 
southeastern Asia, whereas the other taxa distribute in northern China. The Cladiucha species are 
monophagous leaf feeders of Magnolia and Manglietia of Magnoliaceae, a primitive family of 
Angiospermae. Tremex larvae bore into the stems of many species of Salicaceae of Rosopsida and 
Betulaceae and Fagaceae of Hamamelidopsida. The Trichiosoma species feed on the leaves of Salix 
spp. and Populus spp. of Salicaceae, Betula spp. of Betulaceae and Prunus spp. and Sorbus spp. of 
Rosaceae of Rosopsida. The Megalodontes larvae feed on the leaves of Sibiraea spp. of Rosaceae. 
Cladiucha larvae are gregarious, whereas the other aforementioned larvae are all solitary. 

Phylogenetic relationship 

Different substitution patterns were observed among the Symphyta lineages, and all the datasets 
appeared to be biased. However, the AA dataset was the least likely to violate the phylogenetic 
assumptions, and the heterogeneity model performed best in restoring the monophyly of Apocrita and 
Hymenoptera and supporting the formation of the monophyletic group of Apocrita and the Orussidae. 


Thus, we used tree No. 3 as a reference phylogeny (Fig. 9) 


Previous research studies have not considered the heterogeneity of the Symphyta sequence. Some 
of the novel methods that explicitly account for bias can overcome the yield problem caused by 
heterogeneity. However, discrepancies between such inferences within Tenthredinidae still exist. 

The relationships between more heterogeneous taxa, such as Vespina, Siricoidea, Xiphydriidae, 
and Cephidae, have not been resolved in the homogeneity model. However, the heterogeneity model 
recovered Orussidae within Vespina with Apocrita. 

All four schemes of the heterogeneity models support the monophyly of Pamphilioidea. However, 
the relationships between Cephidae remain unresolved. 

Given that all three schemes of the heterogeneity model produced topology No. 2 of 
Tenthredinidae, we speculated that no advantage was obtained in using the heterogeneity model to 
determine the relationship between non-heterogeneous species. 

The ML of all of the datasets did not support the monophyly of Pamphilioidea or perform better 
than BI in restoring the monophyly of Unicalcarida. In determining the internal relationship of 
Cephidae, the problems solved by BI for the same dataset could not be resolved by ML. Thus, we 
speculate that ML is inferior to BI in terms of tree construction, which is consistent with the results of 
many previous studies. 

In contrast to the conclusions drawn from a study of Coleopteral?, we found that the 
mitogenomic data containing RNA reduced the support of the Unicalcarida branches to varying 
degrees. The multiple branches inside the Cephidae also collapsed under RNA addition. We 
speculated that simply using the RNA primary sequence alignment results introduces more noise into 
the tree construction. Exploring the use of the secondary structure to construct phylogeny may lead to 
breakthroughs in the rational use of RNA genes. 

Systematic placement of Cladiucha 

Based on morphological phylogeny, Wei & Nie!!! proposed the Tenthredinoidea s. str. System, 
which should form a framework for subsequent research of Symphytan phylogeny. Some of the results 
were subsequently confirmed by molecular datalé!,, morphological datal?^6, or combinations (8. 
Nevertheless, Cladiucha has not sampled in the cladistic framework except for Wei & Niel! !l. 

Wei & Nie (1998) categorized Cladiucha as neither Allantinae nor Tenthredininae and established 
Megabelesinae, which included Cladiucha, Megabeleses, and two relative genera, based on 
morphological phylogeny. Our mitochondrial genomic data place Cladiucha as a sister to either 
Tenthredininae (topology No. 5) or Fenusinae (topology No. 2 and No. 4). 

The monophyly of Allantinae + Tenthredininae + Megabelesinae (topology No. 5 and No. 6) is 
strongly supported by many apomorphic characters, such as a derived upper head (for example, the 
differentiated frons), a large and broad mesepimeron covering the post-thoracic spiracle, a long and 
multisegmented lancet and a lance with a very short radix. However, the evidence from rRNA went 


against it. The violation may have occured due to the limited taxon sampling. 


Placing Cladiucha as a sister to Allantinae is not supported by our results. In addition to the 
mitochondrial genome data, the following evidence supports this result. The adult morphological 
characters support that Allantinae is a sister group of Tenthredininae; for example, both groups share 
derived mandibles, an enlarged head and a broadened anterior lobe of the pronotum, whereas, in 
Megabelesinae, the mandibles are simple and bidentate, the head is not distinctly enlarged behind the 
eyes, and the anterior lobe of the pronotum is very narrow. Megabelesinae species share a peculiar 
apomorphic character: the female lance has a long and broad membranous lobe and approaches the 
upper 0.3—0.5 of the lancet, which supports that Megabelesinae may have been differentiated and 
isolated from the Allantinae-pattern ancestor over a considerably long period. 

Dated phylogeny of Symphyta 

Li et al.55! depicted the phylogeny of angiosperm, with nearly 3000 chloroplast genomes from 
species representing all 64 orders. Molecular clock dating supports the earliest origin of flowering 
plants, dating the crown group to the Later Triassic period (~209 Ma), which corresponds with the 
timing of the diversification of the largest superfamily in Symphyta, Tenthredinoidea (marked as a 
green block in Fig. 11). Subsequently, the occurrence of the spectacular diversification of core 
angiosperms may arise during the early stages of mesangiosperm evolution (~164 Ma-159 Ma), which 
coincides with the rapid radiation of Tenthredinidae, which is the most diverse lineage in Symphyta. 

All of the known host plants of the Megabelesinae species are plants of Magnoliaceae, a very 
primitive angiosperm family. The magnoliid orders diverged in the Lower Cretaceous (~140-132 Ma), 
as estimated in Li et al./69), or even older!®!, Either the former ог the latter is far earlier than the split of 
Megabelesinae with others. This once again proved that more taxon needs to be sampled in the 
phylogeny. The species differentiation within Cladiucha corresponds with the genus and species 
differentiation within Magnolia!®’! (marked as a green block in Fig. 11). 


5. Conclusions 


1. tRNA gene rearrangements have been found in the mitochondrial genomes of Cladiucha. Two 
species both have the highest A+T content in the Tenthredinidae and, counting only the third codon, in 
the Symphyta. However, neither the ID nor the AliGROOVE showed heterogeneity. 

2. Heterogeneity models can better describe the kinship between heterogeneous species. However, 
there is no obvious substantive improvement in analyzing the relationship of the species that do not 
show the heterogeneity. 

3. The conserved motif in rRNAs could serve as a better potential synapomorphy than the optimized 
alignment. 

4. Tenthredinidae is a monophyletic group. Within this family, ((((Tenthredininae + Fenusinae) + 


Allantinae) + Megabelesinae) + Nematinae) is probably the most acceptable cladogram for the 
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phylogeny of Tenthredinidae, which is also supported by a comparative analysis of the RNA 


secondary structure. 
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Figure Legends 


Fig. 1 Bubble chart of the AT- and GC-skew and A+T% of mitochondrial genomes of Symphyta based on PCG+rRNA data set. (a) 
contains the bubble charts of 44 species in this study, and can be split into (b) to (f) according to different lineages; (b) contains 
Xyelidae (red) and Unicalcarida except for Cephidae (pink); (c) contains Pamphilioidea (orange) and the outer groups (purple); (d) 
contains Tenthredinidae (blue); (e) contains Argidae (light green), Pergidae (light green), Cimbicidae (deep green); (f) contains 
Cephidae (yellow). The letter above the dot represents the abbreviation of the species name. 

Fig. 2 Three-dimensional scatter-plot of the AT- and GC-skew and A+T % of the mitochondrial genomes of Symphyta based on 
PCG+rRNA (A), PCGI2 (B), 9PCG (C) and PCG (D) data sets. PCG + rRNA data set contains 13 protein coding genes and two 
tRNAs; PCGI2 data set contains 13 protein coding genes (excluding the third site); 9PCG data set contains 9 protein coding genes 
(excluding 4 supersaturated protein coding genes); PCG data set contains 13 protein coding genes. 

Fig. 3 Predicted rrnS secondary structure in the Cladiucha mitochondrial genome. The numbering of helices follows Gillespie et al. 
(2006). Roman numerals refer to domain names. Tertiary inter-actions and base triples are connected by continuous lines. C. 
magnoliae as a basemap and base change among Cladiucha species are presented in circles with red (C. magnoliae) and pink (C. 
punctate) color. Logo is generated from 8 species of Tenthredinidae (C.punctate, C. magnoliae, Tenthredo tienmushana, Allantus 
luctifer, Asiemphytus rufocephalus, Analcellicampa xanthosoma, Birmella discoidalisa, Monocellicampa pruni), each Logo graph is 
composed of stacks of letters, with one stack for each position in the consensus RNA secondary structure and the size of the letters 
represents the degree of similarity between bases. 

Fig. 4 Predicted rrnL secondary structure in the Cladiucha mitochondrial genome. The numbering of helices follows Gillespie et al. 
(2006). Roman numerals refer to domain names. Tertiary inter-actions and base triples are connected by continuous lines. C. 
magnoliae as a basemap and base change among Cladiucha species are presented in circles with red (C. magnoliae) and pink (C. 
punctate) color. Logo is generated from 8 species of Tenthredinidae (C. punctate, C. magnoliae, Tenthredo tienmushana, Allantus 
luctifer, Asiemphytus rufocephalus, Analcellicampa xanthosoma, Birmella discoidalisa, Monocellicampa pruni), each Logo graph is 
composed of stacks of letters, with one stack for each position in the consensus RNA secondary structure and the size of the letters 
represents the degree of similarity between bases. 

Fig. 5 Mitochondrial genome organization of Cladiucha with reference to the ancestral type of insect mitochondrial genomes. Genes 
transcribed from the J- and N-strands are shown in green and orange, respectively. The A+ T-rich region is indicated by blue, and 
tRNA genes are labeled by their single-letter amino acid code. 

Fig. 6 Predicted secondary structures of 22 tRNA genes of Cladiucha. Dashes indicate Watson-Crick base pairs, and dots indicate 
G-U base pairing. 

Fig. 7 A+T content and disparity index chart. The chart on the left shows the AT content of 44 species: AT% , AT-1%, AT-2% and 
AT-3%. On the other side of the chart is the AT content of the nucleotides and amino acid (nt123Sum, nt12Sum, ntl Sum, nt2Sum, 
nt3Sum and AA) for 44 species. 

Fig. 8 Heterogeneity analysis of PCGs and PCGRNA datasets. AliGROOVE heat maps of pairwise sequence comparisons for the 
protein coding genes with the nucleotide datasets (PCG123, PCGI23RNA, PCGI2, PCG12RNA,9PCG123 and 9PCG123RNA) and 
the amino acid datasets (АА and 11AA) for the 44 species. The AHGROOVE graph shows the mean similarity scores between 
sequences. AliGROOVE scores range from - 1 (indicating great difference in rates from the remainder of the data set, ie, red coloring 
implies the significant heterogeneity) to +1 (indicating rates match all other comparisons, ie, blue labeling). 

Fig. 9 Phylogenetic hypothesis resulting under MtArt + CAT from Phylobayes analysis of 11AA data set. Branch colors summarize 
support values (see key). Box charts presented on backbone branches represent clade support in each scheme (see “Scheme support"). 
Methods and datasets of phylogenetic tree are summarized in the upper right corner of the figure. The dotted box refers to the 7 
branching models reconstructed with different methods and datasets. 

Fig. 10 Secondary structure drawings of the homology. The conserved structure of rrnS and rrnL is on the right side of each branch. 
The numbers in parentheses represent the number of homologues in the family. The homologous stem-loop structures on the branch 


are framed in the structural diagram. 
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Fig. 11 Dated phylogeny of Symphyta. The axis on the bottom refers to million years and shows the geological time. The blue bars on 
the nodes represent 95% of high posterior density of divergence times obtained from MCMC Tree. The divergence times of each node 
obtained from MCMCTree analysis was written in black on the nodes. The red dots and the times written in red on the nodes indicate 
the divergence time obtained from MEGAX. The green and yellow shaded areas delineate the arising of Tenthredinoidea and 
Cladiucha, respectively, depending on the minimum or maximum ages. The phylogenetic trees given at the bottom are the angiosperm 


tree from Li et al (2019) and Magnolia tree from Nie et al., (2008), respectively. 
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Tables 


Supplementary Table 1 Summary information of symphytan mitochondrial genomes used in phylogenetic analyses 


Species Family Accesion number References 

Ingroup Cladiucha magnoliae Tenthredinidae This study 
Cladiucha punctata Tenthredinidae This study 
Labriocimbex sinica Cimbicidae MH136623 Yan et al, 2019 
Corynis lateralis Cimbicidae KY063728 Dogan and Korkmaz, 2017 
Trichiosoma anthracinum Cimbicidae KT921411 Song et al., 2016 
Megalodontes cephalotes Megalodontesidae МН577058 Niu et al., 2018 
Megalodontes spiraeae Megalodontesidae МН577059 Niu et al., 2018 
Megalodontes quinquecinctus Megalodontesidae MG923502 Tang et al., 2019 
Analcellicampa xanthosoma Tenthredinidae MH992752 Unpublished 
Allantus luctifer Tenthredinidae KJ713152 Wei et al., 2014 
Asiemphytus rufocephalus Tenthredinidae KR703582 Song et al., 2016 
Monocellicampa pruni Tenthredinidae JX566509 Wei et al., 2013 
Tenthredo tienmushana Tenthredinidae KR703581 Song et al., 2015 
Birmella discoidalisa Tenthredinidae MF197548 Unpublished 
Xyela sp. Xyelidae MG923517 Tang et al., 2019 
Xyela curva Xyelidae unpulished unpulished 
Megaxyela euchroma Xyelidae unpulished unpulished 
Xiphydria sp. Xiphydriidae MH422969 Ma et al., 2018 
Tremex columba Siricidae MH422968 Ma et al., 2018 
Pamphilius sp. Pamphiliidae MG923504 Tang et al., 2019 
Chinolyda flagellicornis Pamphiliidae MH577057 Niu et al., 2018 
Orussus occidentalis Orussidae FJ478174 Dowton et al., 2009 
Arge similis Argidae MG923484 Tang et al., 2019 
Arge bella Argidae MF287761 Du et al., 2018 
Calameuta filiformis Cephidae KT260167 Korkmaz et al., 2016 
Calameuta idolon Cephidae KT260168 Korkmaz et al., 2016 
Cephus cinctus Cephidae FJ478173 Dowton et al., 2009 
Cephus pygmeus Cephidae KM377623 Korkmaz et al., 2015 
Cephus sareptanus Cephidae KM377624 Korkmaz et al., 2015 
Characopygus scythicus Cephidae KX907848 Korkmaz et al., 2018 
Hartigia linearis Cephidae KX907843 Korkmaz et al., 2018 
Janus compressus Cephidae KX907844 Korkmaz et al., 2018 
Pachycephus cruentatus Cephidae KX907845 Korkmaz et al., 2018 
Pachycephus smyrnensis Cephidae KX907846 Korkmaz et al., 2018 
Syrista parreyssi Cephidae KX907847 Korkmaz et al., 2018 
Trachelus iudaicus Cephidae KX257357 Korkmaz et al., 2017 
Trachelus tabidus Cephidae KX257358 Korkmaz et al., 2017 
Perga condei Pergidae AY787816 Castro and Dowton, 2005 
Taeniogonalos taihorina Trigonalidae NC027830 Wu et al., 2014 
Parapolybia crocea Vespidae KY679828 Peng et al., 2017 

Outgroup Paroster microsturtensis Dytiscidae MG912997 Hyde et al., 2018 
Neopanorpa phlchra Panorpidae FJ169955 Unpublished 
Neochauliodes parasparsus Corydalidae KX821680 Zhao et al., 2017 


Anopheles gambiae Culicidae L20934 Beard et al., 1993 


Table 1 Summary of mitochondrial genome of Cladiucha magnoliae and C. punctata 


Cladiucha magnoliae 


Cladiucha punctata 


ne ducas Start Stop Length(bp) Start codon Stop codon  Anticodon IGN Start Stop Length(bp) Start codon Stop codon Anticodon IGN 
trnM J 1 69 69 CAT 1 69 69 CAT 

trnQ J 79 147 69 TTG 9 79 147 69 TTG 9 
trn] J 156 222 67 GAT 8 170 236 67 GAT 22 
nad2 J 248 1,300 1,053 ATG TAA 25 287 1,336 1,050 ATG TAA 50 
trnW J 1,306 1,372 67 TCA 5 1,351 1,417 67 TCA 14 
trnC N 1,372 1,443 72 GCA -1 1,417 1,491 75 ССА -1 
trnY N 1,462 1,526 65 GTA 18 1,501 1,566 66 GTA 9 
coxl J 1,535 3,073 1,539 ATT TAA 8 1,573 3,108 1,536 ATT TAA 6 
trnL2 J 3,100 3,168 69 TAA 26 3,154 3222 69 TAA 45 
cox2 J 3,169 3,849 681 ATT TAA 0 3,223 3,903 681 ATT TAA 0 
trnK J 3,854 3,923 70 CTT 4 3,906 3,975 70 CTT 2 
ттр Ј 3,924 3,987 64 СТС 0 3,976 4,039 64 СТС 0 
atps J 3,988 4,149 162 ATT TAA 0 4,040 4201 162 ATT TAA 0 
atp6 J 4,143 4,820 678 ATG TAA -7 4,195 4,872 678 ATG TAA -7 
cox3 J 4,820 5,608 789 ATG TAA -1 4,872 5,657 786 ATG TAA -1 
trnG J 5,614 5,678 65 TCC 5 5,663 5,727 65 TCC 5 
nad3 J 5,679 6,032 354 ATA TAA 0 5,728 6,081 354 ATA TAA 0 
trnA J 6,041 6,106 66 TGC 8 6,099 6,167 69 TGC 17 
trnN J 6,146 6,213 68 GTT 39 6,207 6,275 69 GTT 39 
trnS1 J 6,214 6280 67 GCT 0 6,276 6,342 67 GCT 0 
trnE J 6,288 6,352 65 TTC 7 6,350 6414 65 TTC 7 
ΠΝ Ν 6,351 6,421 71 TCG -2 6,413 6,480 68 TCG -2 
ΠΠ Ν 6,429 6,495 67 GAA 7 6,488 6,555 68 GAA 7 
nad5 N 6,526 8,244 1,719 ATT TAA 30 6,578 8,296 1,719 ATT TAA 22 
trnH N 8,245 8,311 67 GTG 0 8,297 8,364 68 GTG 0 
nad4 N 8,312 9,656 1,345 ATA T 0 8,365 9,709 1,345 ATA T 0 
nad4L N 9,653 9,946 294 ATT TAA d 9,706 9,999 294 ATT TAA -4 
trnT J 9,949 10,019 71 ТСТ 2 10,002 10,069 68 ТСТ 2 
trnP N 10,020 10,086 67 TGG 0 10,070 10,136 67 TGG 0 
ND6 J 10,088 10,600 513 ATA TAA 1 10,138 10,650 513 АТА ТАА 1 
cob J 10,600 11,733 1134 ATG TAA -1 10,650 11783 1,134 ATG TAA -1 
trnS2 J 11,734 11801 68 TGA 0 11,784 11,851 68 ТСА 0 
nadl N 11,840 12,793 954 АТА ТАА 38 11,886 12,839 954 АТА ТАА 34 
trnL1 N 12,794 12,861 68 TAG 0 12,840 12,907 68 TAG 0 
rmL N 12,862 14,205 1,344 0 12,908 14,248 1,341 0 
trnV N 14,206 14,271 66 TAC 0 14,249 14,315 67 ТАС 0 
ттт N 14,402 15,214 813 130 14,467 15,301 835 151 


€^6L000't00202:AIxeuiuo 


